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ABSTRACT 

The dynamical structure of the phase space of the Pluto-Charon system is studied 
in the model of the spatial circular restricted three-body problem by using numerical 
methods. With the newly discovered two small satellites S/2005 PI and S/2005 P2, 
the Pluto-Charon system can be considered as the first known binary system in which 
celestial bodies move in P-type orbits. It is shown that the two satellites are in the 
stable region of the phase space and their origin by capture is unlikely. Also the large 
mass parameter allows to regard the satellites as a model of a new class of exoplanets 
orbiting around stellar binary systems. 

Key words: celestial mechanics ~ planets and satellites: general - Kuiper Belt - 
methods: numerical 



1 INTRODUCTION 

In 1930 C. Tombaugh discovered Pluto, the ninth planet of 
the Solar system. Pluto's fir st moon, Charon was found by 
IChristv fc Harrinertonl lll978f) . The Pluto-Charon system is 
remarkable, since in the Solar system Charon is the largest 
moon relative to its primary, with the highest mass-ratio 
0.130137. Subsequent searches for other satellites around 
Pluto had been unsuccessful until mid M ay 2005, when two 
small satellites were discovered (Stern et al. 2005) provision- 
ally designated as S/2005 PI and S/2005 P2 (henceforth 
PI and P2). With this observation Pluto became the first 
Kuiper Belt object known to have multiple satellites. These 
new satellites are much smaller than Charon, with diame- 
ters 61-167 km (PI) and 46-137 km (P2) depending on the 
albedo. Both satellites appear to be moving in nearly circu- 
lar orbits in the same orbital plane as Charon, with orbital 
periods 38 days (PI) and 25 days (P2). 

From a dynamical point of view the Pluto-Charon sys- 
tem corresponds to such a binary system whose mass pa- 
rameter is approximately one tenth. The phase space of bi- 
naries and the Pluto-Charon system can be studied simul- 
taneously. To survey the phase space of binaries is a fun- 
damental task, since more and more exoplanetary systems 
are being discovered. The great majority of exoplanets have 
been observed axound single stars, but more than 15 planets 
are already known to orbit one of the stellar components in 
binary systems (this type of motion is referred to as satel- 
lite or S-type motion). Other observational facts also favour 
the study of the phase space of binaries. Observations show 
that 60% of t he main sequence stars are in binary or mul- 
tiple systems jPuauennov fc MavorllToQlh . These facts give 
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grounds for the investigation of the stability properties of 
planetary orbits in binaries. 



With the increasing number of extrasolar planetary 
systems around single stars, general stability studies of 
such systems have become of high concern. Many in- 
vestigations were conducted for individual systems, espe- 
cially for those whi ch harbour multiple planets like Ups 
Andromedae see e.g. (|LissauBi| ^99l: [Lauehli n fc Chambers! 



I2OOII: IStepinski. Malhotra fc Black! l2000t). Gliese 876 see 
e.g. (iKinoshita fc N akai' 2001': 'Hadiidcmctriou! !2002!) or 
HD82943 see e.g. iHadiidcmo triou (2002). There are also 
studies on planetary orbits in binaries. The so far discov- 
ered planets in binaries move in S-type orbits. Theoretically 
there is another possible type of motion, the so-called plan- 
etary or P-type, where a planet moves around both stars. 
There are several studies on S- and P-type motions using 
the model of the planar ell i ptic restricted three-body prob- 
lem see e.g. l!Dvorak ''l984'. '1986"; !Holman fc Wiegerl!ll99gl : 
!Pilat-Lohinger fc Dvo rak 2002) and references therein. The 
three-dimension al case, that is the effect of the i nclina tion 
was studied by !Pilat-Lohinger. Funk fc Dvorak! i2003!) for 
P-type orbits in equal-mass binary models. 



The main goal of this paper is to investigate the dy- 
namical structure of the phase space of the Pluto-Charon 
system which can be considered as the first known binary 
system in which celestial bodies, namely PI and P2 move 
in P-type orbits. In Section 2 we describe the investigated 
model and give the initial conditions used in the integra- 
tions. The applied numerical methods are briefly explained 
in Section 3. The results are shown in Section 4. Section 5 
is devoted to some conclusions. 
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2 MODEL AND INITIAL CONDITIONS 

To study the structure of the phase space of the Pluto- 
Charon system we apphed the model of the spatial circular 
restricted three-body problem. We integrated the equations 
of motion by using a Bulirsch-Stoer integrator with adap- 
tive stepsize control. The orbits of the primaries were con- 
sidered circular and their mutual distance A was taken as 
unit distance. The orbital plane of the primaries was used as 
reference plane, in which the line connecting the primaries 
at t = defines a reference x-axis. We assume that the line 
of nodes of the orbital plane of the massless test particle 
(i.e. PI or P2) coincides with the a:-axis at t = 0, thus the 
ascending node = 0° . The pericenter of the test particle's 
orbit is also assumed to be on the a:-axis at t = 0, thus the 
argument of the pericenter w = 0°. Though PI and P2 are in 
the orbital plane of Charon, still we study the problem more 
generally by considering the effect of non-zero inclinations 
on the orbital stability. Thus our results are applicable to a 
wider class of satellite or planetary systems around binaries 
for the mass parameter /i — m2/(mi +m2) = 0.130137, cor- 
responding to the Pluto-Charon system (mi and m2 being 
the mass of Pluto and Charon, respectively). 

To examine the phase space and the stability properties 
of P-type orbits, we varied the initial orbital elements of the 
test particle in the following way (see Table 0: 

• the semimajor axis a is measured from the barycentre 
of Pluto and Charon and it is varied from 0.55 to 5 A (to 4 
A in Fig.El with stepsize Aa = 0.005 A, 

• the eccentricity c is varied from to 0.3 with stepsize 
Ae = 0.05 (Ae = 0.002 in Fig.|3l, 

• the inclination i is varied from 0° to 180° with stepsize 
Ai = 1°, 

• the mean anomaly M is given the values: 0° , 45° , 90° , 
135°, and 180°. 

The above orbital elements refer to a barycentric refer- 
ence frame, where the mass of the barycentre is mi -f- m2. 
By the usual procedure we calculated the barycentric coordi- 
nates and velocities of the test particle and then transformed 
them to a reference frame with Pluto in the origin. In the 
numerical integrations we used the latter coordinates and 
velocities. 

In addition to P-type orbits, S-type orbits were also 
integrated to explore the phase space between Pluto and 
Charon. The initial conditions for S-type orbits are also 
listed in Tabled these are the same both for orbits around 
Pluto and around Charon. 

In total almost six million orbits were integrated for 
10'' Charon's period (hereafter Tc) and approximately 500 
thousand for 10^ Tc- 



3 METHODS 

(i) The relative Lyapunov indicator (RLI). To de- 
termine the dynamical character of orbits we used three 
methods. The method of the relative Lyapunov indicator 
(RLI) was introduced by ISandor. Erdi fc Efthvmiopoulod 
i2000fl for a particular proble m, and its e fficiency was 
demonstrated in a later paper llSandor et alJ 120041 for 2D 
and 4D symplectic mappings and for Hamiltonian systems. 



This method is extremely fast to determine the ordered or 
chaotic nature of orbits. 

The method is based on the idea that two initially nearby 
orbits are integrated simultaneously and the evolution of 
their tangent vectors are also followed. For both orbits the 
Lyapunov characteristic indicator (LCI) is calculated and 
the absolute value of their difference over time is defined as 
the RLI: 

RLl(t) = i|LCl(2;o) - LCl(a::o + Aa::)|, (1) 

where xo is the initial condition of the orbit and Ax is the 
distance of the nearby orbit in the phase space. The value 
of the RLI is characteristically several orders of magnitude 
smaller for regular than for chaotic orbits, thus they can be 
distinguished easily. 

(u) The maximum eccentricity method (MEM). For 
an indication of stability a straightforward check based on 
the eccentricity was used. This action-like variable shows 
the probability of orbital crossing and close encounter of 
two planets, and therefore its value provides information on 
the stability of orbits. We examined the behaviour of the 
eccentricity of the orbit of the test particle along the inte- 
gration, and used its largest value ME as a stability indi- 
cator; in the following we call it the maximum eccentric- 
ity method (hereafter MEM). This is a reliable indicator 
of chaos, since the overlap of two or more resonances induce 
chaos and large excursions in the eccentricity. We know from 
experience that instability comes from a chaotic growth of 
the eccentricity. This simple check was already used in sev- 
eral stability investigations, and was found to be a powerful 
indic ator of the stability character of orbits iDvorak et alJ 
200^ ISiifi. Dvorak fc Freistettei^l2005^ . 

(iii) The maximum difference of the eccentricities 
method (MDEM). We developed this new method and ap- 
plied it for the first time in this investigation. Two initially 
nearby trajectories emanating from a chaotic domain of the 
phase space will diverge according to the strength of chaos. 
The divergence manifests itself in the differences between the 
eccentricities of the orbits and in the angle variables. The 
more chaotic the system is, the faster the difference in the 
eccentricities grows. This difference is sensible to the vari- 
ations around the running average of the eccentricity and 
depends also on the position along the orbit. Thus if the 
positions along the two orbits change chaotically, the eccen- 
tricities of the two orbits also behave differently and their 
momentary differences can be large even if the average value 
of the eccentricity of each orbit remains small. This method 
characterises the stability in the phase space, whereas the 
MEM does it in the space of orbital elements. We define the 
stability indicator MDE as: 

MDE(t) = max|e(t,a::o) - e{t,xo + Ax)\, (2) 

where xo is the initial condition of the orbit and Aa; is the 
distance of the nearby orbit in the phase space. The method 
of the MDE has the advantage with respect to the MEM that 
in the case of chaotic orbits the MDE grows more rapidly 
than the ME, and while the difference between the ME for 
regular and chaotic motions is only 1-2 orders of magnitude, 
this can be 4-7 orders of magnitude for the MDE and there- 
fore can be detected more easily. 
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Table 1 . In the first three rows the orbital elements from unrestricted fits (epoch = 2452600.5) are listed 
C Buie e t al,r20o3): a, e, i, uj, n and M denote the semimajor axis, eccentricity, inclination, argument of the 
pericenter, longitude of the ascending node, and mean anomaly. In the last column the orbital periods are 
given in days. The orbital elements for P- and S-typc orbits are given with the corresponding stepsizes. 



Object 


a [A] 


e 


i[°] 


con 




M[°] 


T [day] 


Charon 
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257.946 


6.387 


S/2005 P2 


2.487 


0.0023 
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267.14 
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4 RESULTS 

We show the results of our investigations in Figs. - 01 
These were obtained as follows. By varying the initial or- 
bital elements as described in Section 2, we performed the 
integration of each orbit for five different initial values of 
the mean anomaly: M = 0°, 45°, 90°, 135°, and 180°. For 
each M the indicators I^'^\a,e,i) were determined, where 
7^^^' stands for RLI, ME, and MDE, respectively. Any value, 
plotted in the figures, is an average over M: 



4:1 5:1 6:1 



I(a,e,i) = - \ /(*^' (a, e, i) . 
5 ^-^ 

M 



(3) 



We note, that this averaging in the case of the RLI and the 
MDE emphasises the chaotic behaviour of an orbit, while in 
the case of the ME is not so drastic. This is due to the fact, 
that the difference in the RLI or in the MDE for regular 
and chaotic orbits is several orders of magnitude, while in 
the case of the ME it is only one or two orders of magnitude. 

The three methods are not equivalent, however they 
complete each other. For example, the ME of the Earth is 
small, indicating stability, although we know from numerical 
experiments that in fact the Earth is moving on a chaotic 
trajectory with a small but nonzero Lyapunov exponent. 
Therefore the ME detects macroscopic instability (which 
may even result in an escape from the system), whereas 
the RLI and the MDE are capable to indicate microscopic 
instability. 

In most of the simulations the j'*^' values were calcu- 
lated for 10^ Tc- To decide whether this time interval is 
enough to map the real structure of the phase space, several 
test runs were done for a much longer time span, for 10^ Tc- 

Fig.Qsummarises the results of the simulations for 10^ 
Tc in the a, i plane, where one can see the structure of the 
phase space: black for chaotic and light grey for ordered 
motion. In these simulations e = was taken for the initial 
eccentricity of the test particle's orbits. 

In general, all the three methods provide the same 
global structure, however, there are also some details that do 
not quite coincide. These differences are only natural, since 
the methods are based on different quantities. In the case 
of the RLI (top panel) the boundary of the chaotic region 
is not so sharp as in the cases of the ME (middle) and the 
MDE (bottom panel). In the RLI map we can see a mixed 
region (grey colour) in the middle, where chaos is already 
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Figure 1. The results of the lO^Tc simulations for e = in the 
a, i plane. Black marks chaotic, light grey means stable motion. 
The white vertical lines indicate the dynamically possible minimal 
pericenter and maximal apocenter distance of P2 (2.41, 2.63 A) 
and PI (3.23, 3.37 A), see Section 4.2 for details. Above the top 
panel the main mean motion resonances are marked. 



present, but manifests itself on a longer time-scale (it would 
need longer integration time to be detected). 

The boundary of the chaotic region slightly extends 
with the increase of the inclination, starting from 2.15 A 
at i = 0°. This is in close agreeme nt with the result of 
iPilat-Lohinger. Funk fc DvorakH2003l) . obtained for /i = 0.5 
(see Fig0 of that paper). It reaches its largest extension at 
i ~ 80°, where it merges with the chaotic region of the 4:1 
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resonance. Increasing the inclination further, the size of the 
chaotic region slowly shrinks. At i ~ 135° the extension of 
the chaotic zone suddenly drops down to a = 1.65 A, then 
it remains nearly constant and its border runs parallel with 
the i-ajcis. Inspecting Fig. the above described structure 
is similar in all panels. 

Several islands appear on the panels, which are con- 
nected with mean motion resonances: 4:1 at Ur — 2.52 A, 
5:1 at Ur = 2.92 A, and 6:1 at = 3.3 A. The first island is 
the most significant and as mentioned above it is connected 
to the large chaotic sea. It is well visible from the figure that 
the higher the order of the resonance, the smaller is the size 
of the corresponding island. We note that the width of the 
islands is larger when the indicators j'*^' are averaged over 
M, than for individual values of M. The widening effect is 
a consequence of a shift of the location of the resonance Ur, 
since slightly depends on M. The location of a resonance 
depends also on the inclination. In the case of the 4:1 res- 
onance flr decreases as i increases, reaching its minimum 
value at about 80°. Further increasing the inclination, a,- 
begins to grow and reaches its initial value at i ~ 160° . This 
dependence of Ur on i is only less visible in the case of the 
5:1 and 6:1 resonances. 

Fig. was obtained from simulations for a time span 
of 10^ Tc- Comparing Fig. [H with the left panel of Fig. H 
(e = 0), obtained from simulations for 10^ Tc, a close agree- 
ment can be observed. The structures described in Fig. 
are already visible in Fig.|5|for the shorter time span. Thus 
we can conclude that the phase space of the Pluto-Charon 
system can be surveyed in a reliable way by using a time 
span of 10^ Tc- 



4.1 The phase space of the Pluto—Charon system 

We investigated the behaviour of P-type orbits systemat- 
ically by changing the initial orbital elements of the test 
particle as described in Section 2 (see also Table 0. Beside 
direct orbits (i < 90°) we studied also retrograd P-type mo- 
tion (i > 90°) of the test particle. All the integrations were 
made for 10^ Tc- The results are shown in Fig. |21 where the 
indicators / are plotted on the a, i plane for different values 
of e. 

In general, the results show an increase of the chaotic 
area for higher eccentricities: for i < 160° the chaotic re- 
gion grows with e. However, the rate of the increase strongly 
varies with the inclination. It can also be seen that the reso- 
nant islands merge with the growing chaotic zone. The most 
striking feature is that the stability of the retrograd P-type 
motion practically does not depend on e. Inspecting Fig. |5| 
it is evident that the border of the chaotic zone for i > 160° 
stay almost constantly at a « 1.7 A. 

Simulations to investigate the phase space between 
Pluto and Charon were also performed for time spans 10^ 
Tc- We examined S-type orbits around Pluto and Charon 
(see the initial conditions in Table 1, where a was measured 
from the respective celestial body). In the case of Pluto 
we found that stable orbits can exist up to 0.5 A, but not 
above this limit. No stable S-type orbits were found around 
Charon. 
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Figure 3. Stability map in the a, e plane. The dark area is un- 
stable, the grey regions are stable. 



4.2 StabiUty of the satelUtes PI and P2 

We have addressed the problem of stability of the recently 
discovered satellites of Pluto. The results are shown in Fig. 
|H] where the values of the MDE (computed for 10^ Tc and 
averaged over for the mean anomalies) are plotted on the 
a, e plane for the planar case {i = 0°). Below a — 2.15 A 
the system is unstable for all e, above a — 2.15 A there is a 
stable region depending on e. The two satellites are situated 
here, in the small rectangles, indicating their dynamically 
possible most probable places of occurance. These rectan- 
gles are defined by the ME, computed in the vicinity of each 
satellite. This means that we took a grid around the present 
value of a of each satellite with a stepsize Aa — 0.005 A, 
Ai — 1.25° in the interval i — — 180°, and with initial 
e = we computed the largest ME,nax during 10^ Tc (in- 
cluding averaging over the five values of M). We obtained 
that MEruax = 0.045 for P2 and MEmax = 0.02 for PI. 
In Fig. 13 these values give the height of the rectangles. We 
computed the possible minimal Vp = a{l — KEmax) pericen- 
ter and maximal Va = a{l + MEmax) apocenter distances of 
the satellites, these are 2.41 and 2.63 A for P2 and 3.23 and 
3.37 A for PI. These values define the horizontal limits of 
the rectangles in Fig. |3 and also the places of the vertical 
lines in Figs. Q and |21 

From Fig. 13 it can be seen that the determined orbital 
elements of the two satellites are well in the stable domain 
of the phase space. If P2 and PI move in the orbital plane 
of Charon, their eccentricities cannot be larger than 0.17 
and 0.31, respectively. The present semimajor axis a of P2 
is very close to the 4:1 resonance, whereas that of PI is 
close to 6:1. The locations of the exact resonances are well 
inside the small rectangles. We presume that these satellites 
probably move in resonant orbits. This could be confirmed 
by new observations. 



5 CONCLUSIONS 

Up to now the Pluto-Charon system is the only known bi- 
nary system which has a relatively large mass-ratio and ce- 
lestial bodies revolve around it in P-type orbits. This circum- 
stance and the high ratio of binary stellar systems among 
stars make it important to study the stability properties 
of P-type orbits in binaries and particularly in the Pluto- 
Charon system. Our investigations show that the stable re- 
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Figure 2. The results of the 10 Tc simulations for e = 0, . . . , 0.3 in the a, i plane. See text for details. 



gion is wider for retrograd than for direct P-type orbits. 
With the increase of the eccentricity the chaotic region be- 
comes larger, and because of it the eccentricities of the two 
satellites, at their present semi-major axis, cannot be higher 
than 0.17 for P2 and 0.31 for PI. Below a = 2.15 A or- 
bits are unstable for all eccentricities, thus no satellite could 
exist here. 

IStern et all (l2005l) has shown that PI and P2 were very 
likely formed together with Charon, due to a collision of 
a large body with Pluto, from material ejected from Pluto 
and/or the Charon progenitor. This is based on the facts 
that PI and P2 move close to Pluto and Charon in nearly 
circular orbits in the same orbital plane as Charon, and they 
are also in or close to higher-order mean motion resonances. 

Our results are also against the capture origin of these 
satellites. Firstly, since the stability region for retrograd or- 
bits is wider, it would have been more probable for the satel- 
lites to be captured into retrograd than for direct orbits. Sec- 
ondly, capture into orbits close to the Pluto-Charon binary 
cannot be with high eccentricity (e > 0.17 and 0.31 at the 
present semimajor axis of P2 and PI), since these orbits be- 
come unstable on a timescale of 10^ Tc- On the other hand, 
for eccentric capture orbits the tidal circularisation time for 
PI and P2 is much l onger than the age of the Solar System 
as lStern et alJ (l2005l) pointed out. 

Taking that PI and P2 were formed together with 
Charon it may be assumed that other satellites were also 
formed from the collision material. However, as Fig.|5|shows, 
to a « 2 j4 binary separation for all inclinations and eccen- 
tricities the orbits become unstable in a short timescale due 
to close encounters with Pluto or Charon. Thus most of the 
ejected material must have been accumulated onto Pluto or 
Charon. How PI and P2 evolved to their current orbits is 
an unsolved problem. 

IStern et~a l. (2005) also suggested that the stochastic 
bombardment of PI and P2 by small Kuiper Belt debris 
can generate transient, dusty ice particle rings around Pluto 
between the orbits of PI and P2. If these particles are co- 
planar with the satellites, their eccentricities must be below 
0.17-0.31 depending on their semimajor axis as Fig.j^shows, 
otherwise they would become unstable in 10^ Tc- 

The existence of S/2005 PI and S/2005 P2 shows that 
there can be a new class of exoplanets which revolve in P- 
type orbits in binary stellar systems. These planetary sys- 
tems can be more similar to our Solar System as compared 
to the known exoplanetary systems, since there cannot be 
large eccentricity orbits in the inner part of these systems. 
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